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Abstract 

The cosmic muon generator CMSCGEN is based on a parametrization of the differential muon flux 
at ground level, as obtained from the air shower simulation program CORSIKA. We present the un- 
derlying ansatz for this parameterization and provide an approximation of the momentum and angular 
distributions in terms of simple polynomials, in the momentum range 3 — 3000 GeV. 



1 Introduction 



The flux <£> of protons hitting the atmosphere is steeply falling with energy, approximately ~ E~ 2 - 7 . This translates 
into a muon momentum spectrum at the surface of the earth falling roughly as p~ 3 . On a flat surface the distribution 
of the zenith angle 8 is to O'^-order proportional to | cos 6\, the exact shape depends on the muon momentum. 

The cosmic muon generator CMSCGEN was adapted from the fast L3CGEN [ 3 1 program which has been 
written for the L3 Cosmics project [ 13 1 some 10 years ago. These programs are based on a parametrization of the 
cosmic muon flux as a function of momentum and zenith angle. Several developments of the last years suggest to 
take a fresh look at this parametrization of the muon flux 0: 

• L3CGEN was successfully used in L3 [6], and since 2006 CMSCGEN has been used[2|. Therefore a pre- 
cise parametrization with known uncertainties is important, now and in future cosmic ray tests. This is in 
particular relevant for collider detectors, for example at the LHC. 

• The air shower program CORSIKA (8), which allows the use of several interaction models, has evolved a 
lot, in particular it is now possible to simulate showers with large zenith angles, while the old version was 
limited to muons with an angle of at least 30° with respect to the horizon. 

• The nuclear interaction models available through CORSIKA have been improved in many ways, and there 
is a brand-new one, EPOS [9|, which is based on recent RHIC data. 

• Precise muon flux measurements [1,6] were made and can be compared to the simulated flux to assess the 
precision of our parametrization. 

We study the muon flux for momenta from 3 GeV to 3000 GcV at the surface of the LHC ring, for all zenith angles. 
This note presents and discusses the parametrization obtained with the CORSIKA version 6.60 (8). 

In the following we first describe the mathematical form of the muon flux parametrization and give the values 
obtained from a CORSIKA simulation for primary protons interacting with the atmosphere, using the default 
interaction models. Then we compare the results to experimental data and to simulations based on other interaction 
models and to the parametrization of 1998 [3|. Finally we discuss the ratio of the flux for positive and negative 
muons and the influence of heavier primary nuclei on the momentum spectrum of the muons. 



2 Muon Flux Parametrization 

We follow closely the procedure outlined in Q : The total muon flux, integrated over the full range of zenith angles 
(vertical to horizontal) and all azimuthal angles, at an altitude of about 500 m is parametrized in the form 

d$ 1 

— ; — = C norm ■ r • s(L) . (1) 
dp p° 

$ denotes the number of muons per time and area hitting a (horizontal) surface area, p is the muon momentum and 
L = log 10 (p/GcV). The flux normalization constant C norm will be discussed later. With this parametrization of 
the steeply falling muon momentum spectrum the function s(L) (which we need to determine) varies rather little 
with L. The definition (Q]) implies that s(L) can be written as 

s{L)~p>—. (2) 

This distribution can be obtained from the histogram of the logarithm L, by applying a weight of (p/Gc-V) 2 = 10 2L 
to each muon generated by CORSIKA. Since the cosmic muon flux falls roughly with p~ 3 , the expression s(L) is 
a slowly varying function of L, and we can approximate it by a polynomial: 

s(L) = a + ai L + ... + a 6 L e (3) 

It turns out that a polynomial of degree six is sufficient for our purposes. 

The distribution of the azimuthal angle tfi is assumed to be flat, the dependence on the momentum dependent zenith 
angle (here: 9 = 180° for a vertically downgoing muon) is parametrized by the normalized function z(c, L) with 
c = cos 6: 

^ ~ z(c, L) = b (L) + h(L) c + b 2 (L) c 2 (4) 
dc 
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We use spherical coordinates with the z axis pointing vertically upwards; thus vertical muons correspond to c = 
— 1, horizontal ones to c = 0; these numbers define the range of c values to be considered. Note that the three 
coefficients bi are not independent from each other since the integral J z(c, L) is one by definition. The coefficients 
bi depend on the muon momentum — we will discuss this later together with the Monte Carlo data. 

The global constant C norm can be derived by comparing the predicted muon flux with the measured one, for 
vertical muons at the reference momentum of 100 GeV, see below. 

Thus the full parametrization has the following simple form: 



dp dc dtp 



4t ■ s(L) ■ z(c,L) 



2tt 



(5) 



3 CORSIKA Simulation 

The primary particles impinging on the atmosphere are assumed to be protons with an energy distribution E~ 2 ''. 
We use CORSIKA d with the model EPOS (version 1.61) JD for high enery interactions (lab energy of 80 GeV 
or more) and GHEISHA (version 2002d) lfl2ll for low energy hadronic interactions. The standard atmosphere is 
used, zenith angles 8 with | cos#| > 0.087 are generated (corresponding to showers 5° above the horizon). In the 
following we limit the cos 8 range to the interval [—1, —0.1]. 

We have generated five samples of 1 million showers each, for different primary energy ranges, covering in total 
all proton energies from 5 GeV to 10 7 GeV. Lower or higher values hardly contribute to the muon momenta we 
are investigating here (3J. 




L = log io (p/GeV) 

Figure 1: Momentum spectrum expressed through function s(L), integrated over zenith angle and azimuth. The 
CORSIKA points are obtained with interaction packages EPOS and GHEISHA. 

Figure Q] shows the the function s(L) as simulated for the total muon flux, together with the fitted polynomial 

s(L) = -1 + 6.2218 • L - 13.940 • L 2 + 18.164 • L 3 ~ 9.2278 • L 4 + 1.9923 • L 5 - 0.15643 • L 6 (6) 

Note that we can normalize the coefficients of s(L) arbitrarily; here we have done it such that the modulus of 
the first coefficient is one. The fit reproduces the simulated points with an accuracy of the order of 5%, in the 
momentum range from p — 3 GeV (L w 0.5) to p — 3000 GeV (L w 3.5). At the lowest L-bins the fit degrades 
somewhat (hardly visible in FigureQ] but the deviations reach up to 20%) - we will come back to this later. 

The cos 8 distribution z(c) is shown for two momentum ranges (near 10 GeV and around 1000 GeV) in Figure|2] 
for —1 < c < —0.1. Again the CORSIKA points are approximated quite accurately by the parametrization. It is 
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Figure 2: Distribution z(c) with c = cos 6 (zenith angle) for two different muon momenta. The CORSIKA points 
are obtained using a combination of the interaction packages EPOS and GHEISHA. 

evident from Figure [2] that the momentum dependence of the zenith angle distribution must be taken into account. 
Note that for an isotropic muon flux we would expect z ~ | cos 0\, due to the flat surface we consider here (and not 
a spherical one). 

The L dependence of the coefficients in © can again be parametrized by a polynomial. We obtain: 

bo(L) = 0.6639 - 0.9587 • L + 0.2772 • L 2 
b x (L) = 5.820 — 6.864 • L + 1.367 • L 2 

b 2 (L) = 10.39 — 8.593 • L + 1.547 • L 2 (7) 



Technical remarks: please note that with these fitted curves the normalization requirement 



-0.1 



z(c, L) dc = 1 



(8) 



is respected only approximately; consequently the as given above must be renormalized as a function of L before 
z(c, L) can be calculated. Also, for near horizontal showers and small momenta the flux becomes very small and 
the parametrization can result in negative values; this region in (L, c) should be cut away. 



4 Comparison with Experimental Results 

The flux measurements of the last century are summarized in the compilation 1 1 1. The most precise results have 
been obtained by the L3 collaboration in 2004. 

In Figure[3] we compare the measured vertical flu)0 to our parametrization of s v (L), which we have obtained in 
a similar way as s(L) (equation [6]) but selecting only (near-)vertical muonQ We show the ratio instead of the 
flux itself in order to be more sensitive to potential deviations in the distributions. Note: since we have not yet 
introduced an absolute normalization, we have arbitrarily set the ratio to 1 near 100 GeV (L = 2). Since the 
compilation [ 1 1 refers to the flux at altitude (sea level), we applied the small correction as given in formula (1) in 
to extrapolate to the altitude of 500 m. Since the L3 flux is measured at a similar altitude, 470 m, no correction 
is necessary. The agreement between parametrization and measurements is remarkably good — only at momenta 

^ only the vertical flux has been measured by several detectors 

2 ' s v (L) looks similar to s(L) as shown in FigureQ] but the peak is shifted slightly to the left. 
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Figure 3: Ratio of measured and parametrized vertical muon fluxes. The small open circles indicate the uncertain- 
ties of the curve with big open circles, obtained from the compilation [ 1 1. The black circles show the results of the 
L3 measurements (6J. 



above 500 GeV deviations become visible. 

The measured and parametrized cos(6>) distributions are compared in Figure [4] Note that the L3 points have 
been normalized 'by eye' to the parametrization. The overall agreement is good, just for (near-)vertical muons 
(c = —1 ... — 0.95) the line is a couple of percent too high for low momenta. 



5 Comparing Predictions of Different Models 

Finally we compare our parametrization as obtained with EPOS and GHEISHA to the predictions by other models. 

First we exchange GHEISHA against the other low energy hadronic interaction model available in CORSIKA, 
FLUKA (version 2006.3b) iflOl , and we keep EPOS. In a second simulation run we keep GHEISHA, but exchange 
EPOS against QGSJET (version 11-03) ATI . The result is shown in Figure [5] in terms of the ratio of the fluxes 
(zenith angle integrated) predicted by the different models. Again, we have arbitrarily normalized the ratio to 1 
near 100 GeV (L — 2). Differences are clearly visible. The two alternative models introduced here differ from 
each other by up to 30% in the momentum range 3 — 3000 GeV. Thus the models are not yet good enough to predict 
the flux with the same precision as obtained in the recent measurements. Our reference model, EPOS+GHEISHA, 
lies in between the other two models. Since the experimental data are in agreement with the reference model for 
momenta between 10 and 500 GeV, these differences are only relevant at low and high momenta. The discrepan- 
cies here indicate the uncertainties intrinsic in these models. 

How does our new parametrization in equation (O compare to the old one from 1998 [3| ? Figure [6] shows the 
ratio of the fluxes, integrated over | cos6*| > 0.4, the zenith angle range used in the old parametrization. The 
differences are quite large, the new momentum spectrum is significantly harder. Actually, looking at Figures[5]and 
|6l there seems to be a general trend: the newer models (EPOS and FLUKA) give more muons with relatively high 
momenta . . . Clearly, the discrepancies between the old and new parametrization make an update of CMSCGEN 
mandatory. 

Figure [7] compares the parametrizations for the cos 8 distribution. Shown is the angular range and the momentum 
values as displayed in the original note [3|. The agreement is satisfactory; the small difference in shape can be 
attributed to the fact that before a linear approximation was used, while we now use a quadratic polynomial. 
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Figure 4: Distribution z(c) as measured by L3 for two different momentum values, in comparison with the 
parametrization. The 1600 GeV graph is multiplied by 1/2 for better visibility. 
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Figure 6: Ratio of fluxes for old and new parametrization. 
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Figure 7: cos 6 distribution — comparison of new and old parametrization 
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6 Muon Flux Normalization 

For muons of momentum p = 100 GeV the vertical flux has been measured to [ 1 1 



and © 



dp d cos 8 < 



d$ 
dp d cos t 



(2.59 ± 0.18) ■ 10~ 3 irT 2 s" 1 GeV -1 sr" 1 



= (2.63 ± 0.06) • 10~ 3 m" 2 s" 1 GeV -1 sr" 1 



(9) 



(10) 



Combining the two results leads to a flux of vertical muons with the following value, dominated by the L3 mea- 
surement: 



dp dcos 8 ■ 



(2.63 ± 0.06) ■ 10~ 3 m~ z s" 1 GeV sr 



(11) 



To calculate C norm we first compute the product 

P=l. s (L)-z(c,L)-^- (12) 
p a Zir 

for L = 2 and c = — 1 using our parametrizations © and ©: 

P = 10~ 6 GeV~ 3 • 7.094 • [-0.1653 - 2.7882 • (-1) - 0.6948 • (-1) 2 ] • — = 2.177 • 10~ 6 GeV -3 (13) 

In order to compensate for the small discrepancy at c = —1 between L3 data and the parametrization of z(c), 
which occurs at 46 GeV, see Figure |2] and Figure @] but also — to a lesser extent — at 100 GeV, we introduce a 
'fudge factor' of 1.05 to optimize the overall agreement of the parametrization with the L3 results: 



P -> P/1.05 = 2.073 • 10" 6 GeV" 



Comparing with <(5j and (fTTT i gives: 



C norm = 1.27 • 10 3 GeV 2 m" 2 s" 1 sr' 



(14) 



(15) 



7 Charge Ratio 

The measurements yield a ratio R of positive and negative muon fluxes of [1 1 

R = 1.268 ±0.028 



and 



resulting in a combined value of 



R = 1.285 ±0.019 



R = 1.280 ±0.016 



(16) 



(17) 



(18) 



for vertical incidence and momenta around 100 GeV. The data are consistent with a charge ratio R which is 
independent of momentum (at least from 10 to 500 GeV) and also of zenith angle, if one stays away from near 
horizontal muons (6] [l4ll . 

The CORSIKA simulations (EPOS+GHEISHA) yield R values varying by ±5% around 1.40 for momenta in the 
range 10 to 1000 GeV, in disagreement with the measurements. 

For the generator CMSCGEN the experimental result of R = 1.28 must be used, i.e. the total muon flux must be 
split into positive and negative muons with the relative fractions of 



/+=66.1% f- 



43.9% 



(19) 
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8 Heavy Nuclei 

It was already shown in (3j that primary helium or iron nuclei yield muon spectra that are similar to those from 
protons. Since anyway the parametrization presented here describes the experimental data well, there was no need 
to repeat such a study here. 

9 Discussion of Uncertainties 

Our parametrizations for z(c) can be used from c = — 1 (vertically downward going muons) to c = —0.1. For 
higher values (nearly horizontal showers) the Jura mountains introduce a <fi dependent absorption of muons, since 
they stand up about 10° above the horizon, seen e.g. from the CMS[7| site. 

For momenta below 3 GeV (L « 0.5) and above 3000 GeV (L ks 3.5) our parametrization is not valid and should 
therefore not be used. Reasons are the instability of the fit for low and high momenta (see Figures[5]and[6]l and the 
lack of reliable experimental data in these regimes (see Figure[3]). 

In the 'central' momentum range from 10 to 500 GeV our parametrization seems to work quite well: the measured 
cosmic fluxes are reproduced at the 5 — 10% level. So we assign an uncertainty (68% confidence level) to the 
absolute differential flux as given by our parametrization of ±7%. 

In the 'fringe' momentum regions (below 10 GeV and above 500 GeV) the uncertainty increases rapidly. Due to 
the missing experimental support at high momenta we assign a relative uncertainty to the flux increasing up to 
±50% at 3000 GeV. This estimate is obtained by extrapolating the uncertainties displayed in Figure [3] At low 
momenta geomagnetic effects become important [ 14 1 and atmospheric and solar influences make the muon flux 
vary with time [1], there is no comparison to experiment either, and the model predictions disagree with each 
other, see above. Therefore we attribute to the low momentum regime a relative uncertainty increasing from ±7% 
at 10 GeV up to ±25% for p = 3 GeV. 
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